Stability of a BEC with Higher-order Interactions near a Feshbach Resonance 
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We consider the stability of an ultracold trapped Bose-Einstein condensate near a Feshbach res- 
onance. Using a modified Gross-Pitaevskii equation that includes higher-order terms and a multi- 
channel model of Feshbach resonances, we find regions where stability can be enhanced or suppressed 
around experimentally measured resonances. We suggest a number of ways to probe the stability 
diagram. Using scattering length zero-crossings huge deviations are founds for the critical particle 
number. Effects are enhanced for narrow resonances or tighter traps. Macroscopic tunneling of the 
condensate is another possible probe for higher-order interactions, however, to see this requires very 
narrow resonances or very small particle numbers. 
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Introduction The stability of Bose-Einstein conden- 
sates (BECs) in ultracold alkali gases is determined by 
the sign of the scattering length a For a < 0, one 
has effectively attractive interactions and the condensate 
will collapse to a dense state when the number of con- 
densed particles N is larger than a critical number N c 
Pi S S S] ■ This has been beautifully demonstrated in 
experiments with 7 Li IB , 87 Rb @, 0] , and recently with 
a dipolar 52 Cr BEC 8]. The findings indicate that the 
theory based on the surprisingly simple Gross-Pitaevskii 
(GP) equation can reproduce and describe most features 
of the experiments. 

The GP equation includes two-body terms through a 
contact interaction which is parametrized by a. This is 
equivalent to a Born approximation but with an effective 
coupling that is obtained by replacing atom by the phys- 
ical scattering length a. However, higher-order terms in 
the expansion of the phase shifts at low momenta, deter- 
mined by the effective range r e , the shape parameter etc., 
give corrections to the simple GP equation. In this paper 
we explore the influence of the effective range term on the 
properties of a BEC. In particular, we show that the crit- 
ical number of condensed atoms depends strongly on the 
higher-order scattering term when the scattering length 
approaches zero (zero-crossing). We also show how the 
macroscopic quantum tunneling (MQT) rate, in which 
the entire BEC tunnels as a coherent entity, can be mod- 
ified for small condensate samples. 

The considered effects depend on a combination of a 
and r e which yield different behavior for wide and nar- 
row Feshbach resonances. Recent measurements on 39 K 
found many both wide and narrow resonances which al- 
low for tuning of a over many orders of magnitude 0, [l(| • 
We therefore consider a selected example from 39 K in or- 
der to elucidate the general behavior for realistic experi- 
mental parameters. 

The paper is organized as follows: We first introduce 
the modified GP equation with effective range depen- 
dence. Using both a variational and numerical approach 



we find a phase diagram describing the stability of the 
BEC. We then consider a Feshbach resonance model in- 
cluding effective range variations. The behavior of the 
critical particle number near a scattering length zero- 
crossing is derived. We discuss MQT and show numeri- 
cally how the rate is modified. We finally discuss other 
possibilities for probing the higher-order interactions. 

Modified GP Equation We assume that the conden- 
sate can be described by the GP equation and we focus 
on the a < attractively interacting case. Since we are 
interested in the ultracold regime, where the tempera- 
ture is much smaller than the critical temperature for 
condensation, we adopt the T = formalism. In order 
to include higher-order effects in the two-body scattering 
dynamics, we use the modified GP equation derived in 
[ll| for which the equivalent energy functional is 
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where m is the atomic mass, V is the external trap, Uq — 
4irh 2 a/m, and <?2 = a 2 /3 — ar e /2 with a and r e being the 
s-wave scattering length and effective range, respectively. 

We are interested in the stability properties of the 
ground-state and we therefore perform a variational cal- 
culation using the mean-field trial wavefunction 
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where q is the dimensionless variational parameter and 
b = \J h/muj is the trap length. The normalization 
is N = J dr|\f(r)| 2 . For simplicity we only consider 
isotropic traps with V(r) = \muJ 1 r 1 . However, the ef- 
fects found should hold for deformed traps as well (along 
the lines of the analysis in [l^] ) ■ The variational energy 
is 
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In Fig.[T]we plot E(q) for different parameters. As shown, 
there are many possibilities for 172 7^ 0, including stable, 
unstable, and metastable systems. We see that the 32 
term modifies the barrier for N\a\/b = 0.5, implying that 
tunneling rates will be altered. For iV|o|/6 = 0.7, the 
#2 = case has no barrier at all, and here the 32 term 
can in fact produce a small barrier on its own. The g~ 5 
dependence of the term means that the effect is small. 
However, the plot clearly shows that a new stability anal- 
ysis is needed. In addition to the variational approach 
we have numerically solved the full time-independent GP 
equation corresponding to eq. (|TJ) . 
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FIG. 1: Energy of a BEC with fixed N\a\/b as function of 
the variational parameter q, i.e. the size of the condensate. 
The higher-order interaction term </2 modifies the height and 
shape of the barrier. 

Phase Stability Diagram To determine the ground- 
state stability one looks for the vanishing of the bar- 
rier towards q = 0. For g 2 = eq. ^ leads to 
N c \a\/b w 0.671 Q. The full integration of the GP equa- 
tion gives N c \a\/b w k , k = 0.5746 HQ. Th ese values 
are indicated by filled points in Fig. O 

In the general (72 7^ case we first take the variational 
energy eq. ([3]) and solve for multiple roots of dE(q)/dq. 
The resulting "phase diagram" is plotted in Fig. [5] (solid 
line). In the upper left region (S) we have complete sta- 
bility of the condensate; only one minimum exists at large 
q ~ 0.8 (see Fig.[T]) and the potential goes to plus infinity 
at small q. In the metastable region (M) a barrier and a 
minimum exists for large o; for gi > another minimum 
exists at small q < 0.2, while for gi < the potential 
goes to minus infinity. In the unstable region (U) the 
barrier has vanished; the potential either has a minimum 
at small q (52 > 0) or no minimum at all (32 < 0). In 
the variational approach the stable and unstable regions 
are connected via the upper right part of Fig. [21 Go- 
ing from (S) to (U) corresponds to an adiabatic change, 
where the macroscopic (q ~ 1) minimum is transferred 
to a microscopic (g <C 1) high-density minimum. 

Next, we numerically solve the full the GP equa- 
tion. The stationary (white) and non-stationary (shaded) 
regions are shown in Fig. [2l Since both stable and 
metastable solutions are considered stationary, the white 
regions cover both (S) and (M). Both our variational and 
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FIG. 2: (color online) Phase stability diagram of a BEC 
with higher-order interactions. The solid line obtained 
from the variational ansatz, eq. ((2j), divides the stable (S), 
metastable (M), and unstable (U) regions. The white re- 
gion indicates where stationary solutions exists for the full 
GP equation. Filled points show the well-known (72 = re- 
sults. Arrows correspond to the values in Fig. [T] The dashed 
line is the fit in eq. ((5}. Two specific 39 K values are chosen 
for zero-crossing (cross) and MQT rates (open square) calcu- 
lations (see text). 



numerical results agrees with the known 32 = results, 
and have similar behavior for small and negative g2/b 2 . 

Feshbach Resonance Model In order to predict effects 
of the higher-order term, we need a realistic model for 
g2 ■ Since g2/b 2 is the relevant parameter, and the trap 
length, b, is usually orders of magnitude larger than 
atomic scales, it is necessary to look for divergences of 
gi- Since 172 depends on a and r e , a Feshbach resonance 
is the obvious mean. The standard single-channel mod- 
els are inadequate since only a is considered. We there- 
fore use a multi-channel model [l3| , which describes both 
a and r e as a function of resonance position Bq, width 
AB, magnetic moment difference between the channels 
A/z, and the background scattering length a^g. In [l3j |. 
only the regime \B — Bq\ <C AB was of interest, yielding 
r e independent of magnetic field, B. In the current work 
we generalize to keep the .B-dependence of both a and r e . 
This gives a = ab g (l — AB/(B — B )) and r e — r e0 /(l — 
(B-B )/AB) 2 , where r e0 = -2h 2 /ma bg AfiAB < 0. We 
thus have a field-dependent 52 (-B) for given values of Bq, 
AB, A/i, and a& 9 . Notice that r e = r e o(l — ab g /a) 2 and 
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Hence 172 diverges when a — > (referred to as zero- 
crossing) or a — > 00 (on resonance). These are well- 
known features for model potentials like the square-well 
and van der Waals interactions. 

As a concrete example we use the extremely narrow 
Feshbach resonance in 39 K found at Bq — 825G, with 



3 




a/(10 2 a ) 
r e /(10 5 a ) 



-60 -40 -20 20 

(B - B )/mG 



40 



GO 



FIG. 3: (color online) Scattering length a, effective range r e , 
and coupling constant g 2 as function of _B-field for the narrow 
39 K Feshbach resonance at Bo = 825 G. The inset shows 
the MQT rate F 2 (and T for g 2 = 0) for N = 242 and 
b = 1.84/im. The cross and open square are as in Fig. [2] 



AB = -0.032G, Afi = 3.92^ B , and a bg = -36a [9|. 
For this resonance we find a large r e o = —2.93 x 10 4 ao. 
The variation of a, r e and gi as function of B is shown 
in Fig. El We use trap length b = 1.84/xm = 3.48 x 10 4 a 
in all calculations unless indicated otherwise. 

Critical Particle Number Near Zero-crossings As 
noted above, large values of 172 are possible at zero- 
crossings (a = 0) since r e diverges here, see e.g. Fig. [3] 
around B — Bq ~ — 32mG. The scattering length is cur- 
rently being tuned with extreme accuracy near such zero- 
crossings in 39 K around 350G with control down to the 
0.06ao level However, this is for a broad resonance 
and the same level of tuning around narrow ones is likely 
much more challenging, particularly in high magnetic 
fields. Another narrow resonance at Bo = 25.85G with 
AB = 0.47G exists in 39 K 0. It will show the same fea- 
tures as the one considered here and the lower magnetic 
field could make tuning easier. 

The advantage is that many particles N c cx l/\a\ can 
be accommodated in the condensate. However, for a — > 
we have 0172 — ► — ?~eo(a&g) 2 /2, i.e. a finite limit. Remem- 
bering that r e0 and V 2 |^| 2 are negative, the last term 
in the energy functional cq. ([I]) also becomes negative. 
Thus larger densities or density fluctuations gives lower 
total energy. This implies less stability and smaller N c 
near a = 0. 

To calculate quantitative effects on A c we focus on the 
critical line between the unstable and metastable regions 
in Fig. [2] for 52 < 0. We fit the dependency as 
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where ki = 14.5, as shown in Fig. [2j With 32(0) given 
by eq. ((!]), A c becomes a function of \a\/b for fixed values 
of abg/b and r e o/b. In Fig. [4] we plot 1/A C as function 
of \a\/b for the 39 K resonance with various trap lengths. 
The curve deviates from the linear 52 — result both 
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FIG. 4: (color online) Critical particle number N c as function 
of scattering length a. The linear dashed (red) line is the 
well-known g 2 = result N\a\/b = fco. The solid (black) lines 
are for a narrow 39 K resonance and trap length b = 1.84/im, 
(bottom), 0.5& (middle), and 0.26 (top). The cross and open 
square are as in Fig. [2] 



at \a\ ~ and \a\ ~ b. Squeezing the trap makes the 
effect even more pronounced. In the limit \a\/b — > the 
critical particle number approaches a finite value, A c — > 
2k b 3 /(ki\r e r l \al g ). Large effects also occur near \a\ ~ b, 
however, here the allowed particle number is very small. 

To exemplify with the 39 K resonance we use B — Bq = 
— 32.3mG (cross in Fig. [3]) giving a = — 0.334ao and 
g-2/b 2 = —0.0461. The critical particle number from 
eq. |5]) is A c = 3.58 x 10 4 . The values are indicated 
by crosses in Fig. [2] and Fig. |U The g2 — prediction is 
A c = 5.975 x 10 4 . In conclusion, we predict 40% reduc- 
tion in A c for macroscopic particle number A c ~ 10 4 . 

Macroscopic Quantum Tunneling Close to a = 00, 
g-2/b 2 is still a small term in typical traps. As the MQT 
rate [l5| is exponentially dependent on the integral under 
the barrier, we expect changes due to non-zero 32 to be 
amplified. A similar point was raised in [12j where the 
MQT rate was shown to increase dramatically close to 
Ac. The MQT rate, T2, can be obtained from field theory 
[171 ] by considering the bounce solution of the effective ac- 
tion in variational g-space with potential E(q)/Ntku. All 
MQT calculations are done within the variational ap- 
proach, eq. ([3]). Since we assume spherical symmetry, we 
use cj = Q = (uJxUjyUJz) 1 / 3, when relating to experiment. 
For comparison, we denote by Tq the rate with g2 = 0. 

We are interested in metastable states with small bar- 
riers and large MQT rates, thus we have to work close to 
A c Below we find that g2/b 2 is only important for 
the rate when A is of order 10 2 or less. Since the three- 
body loss depends on the density to the third power [l|| , 
we expect it to be small in the outer metastable minimum 
and to be large in the inner one. Thus, the physical pic- 
ture is that of a metastable BEC in the outer minimum 
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that knocks on the barrier as a coherent state with a com- 
mon quantum tunneling probability. When it tunnels to 
the inner minimum it swiftly decays as the large increase 
in density amplifies the three-body loss. Thermal fluctu- 
ations can of course also induce MQT and we therefore 
need to operate at very low temperature. An estimate 
of the thermal tunneling rate is given in 17J , and leads 
to the following criterion for thermal fluctuations to be 
suppressed: E{q m ) — E(qo) ^> ksT, where q m and qo 
are the positions of the barrier maximum and the outer 
minimum, respectively. 

For the 39 K resonance we now consider B—Bq = 20mG 
(see Fig. EJ) where g 2 /b 2 = -4.3 x 10" 4 and a = -93.6a . 
The critical particle number becomes N c ~ 249 (see 
Fig. H} . We now pick N = 242 < N c to obtain metasta- 
bility (see Fig. The rates are found to be T 2 = 0.49 
s _1 and Tq = 0.05 s _1 (shown by the open squares in the 
inset of Fig. [3]), i.e. a tenfold enhancement of the MQT 
rate. The temperature must be below 8 nK for thermal 
fluctuations to be small. It is possible to obtain larger 
N with the cost of a smaller effect on the MQT rate. 
However, going to N ~ 10 3 does not seem possible. 

We propose to start from the a > side where the con- 
densate is stable. Preparing a sample with low enough N 
and temperature of a few nK for the effect to be observ- 
able is the very difficult challenge for the MQT scenario. 
A sudden ramp of the magnetic field to the appropriate 
a < is then performed. The system has to retain a 
small barrier so that the MQT rate is considerable and 
to avoid collapse. Then one monitors N as a function of 
time, as was done in e.g. [?|]- The low ./V is also challeng- 
ing for optical imaging of the atomic cloud, but should 
be possible with current techniques. We expect not only 
to see modified MQT rates, but also changes in the decay 
after the condensate tunnels to a high density state; The 
three-body recombination is very sensitive to the density, 
and the inner barrier caused by the g 2 term should affect 
this process as well. 

Other Signatures From the structure of the 172 term 
it is clear that density variations are needed to see ef- 
fects. As we have demonstrated already, tighter traps 
amplify the contribution, and we expect optical lattices 
to do the same. Solitons and vortices are other features 
with density variation. For a simple one vortex state, the 
|\f r | 2 V 2 |'I , | 2 term will leave the core and the asymptotic 
regions unchanged but change the profile in between. 
Work is in progress to extract these changes fljj ]. 

Conclusion We have explored the effect of higher- 
order terms in the Gross-Pitaevskii description of a Bosc- 
Einstein condensate, particularly the effective range cor- 
rection near a Feshbach resonance. Using both a vari- 
ational and numerical approach we find an interesting 
new phase diagram for the stability of a condensate with 
negative scattering length. 

The critical particle number is strongly affected near 
zero-crossings of the scattering length for narrow reso- 



nances. Effects of 40% are found for particle number of 
order 10 4 . These deviations increase for more narrow res- 
onances or tighter traps. We find that the critical number 
will be reduced as the zero-crossing is approached from 
either side since the higher-order term is always attrac- 
tive. The dipole-dipole interaction is usually important 
near zero-crossings but we estimate it to be smaller than 
the higher-order term for narrow resonances and tight 
traps. Macroscopic quantum tunneling is also modified 
by these higher-order corrections and we have discussed 
some experimental conditions for exploring the physics. 
Narrow Feshbach resonances are the best way to isolate 
the effect of tunneling. However, one needs small samples 
of order 10 2 particles and temperatures around 10 nK. 
Other possible experimental signatures are tight traps, 
optical lattices, solitons, and vortices in rotating BECs. 

In a more general sense, the fact that higher-order in- 
teractions can become large and dominant in determining 
the stability properties means that it is no longer a cor- 
rection and that still higher terms might become impor- 
tant. This is the subject of future work. For the moment 
we have demonstrated that experiments targeting the re- 
gions discussed above could find interesting new stability 
properties beyond those that were already understood 
about a decade ago. 

Acknowledgments Discussions and encouragement 
from Aksel S. Jensen is gratefully appreciated. NTZ ac- 
knowledges useful comments from Doron Gazit, as well 
as support from the Villum Kann Rasmussen foundation. 



[1 

[2 

[3. 

[4] 
[5] 

ir 

[8 
[9 
[10 

[11 

[12 

[13 

[14; 
[15 

[16 

[17. 
[18 



F. Dalfovo, S. Giorgini, L.P. Pitaevskii, and S. Stringari, 
Rev. Mod. Phys. 71, 463 (1999). 

P.A. Ruprecht, M.J. Holland, K. Burnett, and M. Ed- 
wards, Phys. Rev. A 51, 4704 (1995). 
A. Gammal, T. Frederico, and L. Tomio, Phys. Rev. A 
64, 055602 (2001). 

G. Baym and C.J. Pethick, Phys. Rev. Lett. 76, 6 (1996). 
C.C. Bradley, C.A. Sackett, and R.G. Hulet, Phys. Rev. 
Lett. 78, 985 (1997). 

J. L. Roberts et al, Phys. Rev. Lett. 86, 4211 (2001). 

E.A. Donley, Nature 412, 295 (2001). 

T. Lahaye et al, Phys. Rev. Lett. 101, 080401 (2008). 

C. D'Errico et al, New J. Phys. 9, 223 (2007). 

G. Roati et al., Phys. Rev. Lett. 99, 010403 (2007). 

A. Collin, P. Massignan, and C.J. Pethick, Phys. Rev. A 
75, 013615 (2007). 

M. Ueda and A.J. Leggett, Phys. Rev. Lett. 80, 1576 
(1998). 

CM. Bruun, AD. Jackson, and E.E. Kolomeitsev, Phys. 
Rev. A 71, 052713 (2005). 

M. Fattori et al, Phys. Rev. Lett. 101, 190405 (2008). 
A nice introduction to MQT with references is found in 

D. Mozyrsky et al, Phys. Rev. A 76, 051601(R) (2007). 
S.K. Adhikari, Phys. Rev. A 66, 013611 (2002). 

H. T.C. Stoof, J. Stat. Phys. 87, 1353 (1997). 
M. Th0gersen and N.T. Zinner, in preparation. 



